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The detection of signals with varying frequency is important in many areas of physics and 
astrophysics. The current work was motivated by a desire to detect gravitational waves from the 
binary inspiral of neutron stars and black holes, a topic of significant interest for the new generation 
of interferometric gravitational wave detectors such as LIGO. However, this work has significant 
generality beyond gravitational wave signal detection. 

We define a Fast Chirp Transform (FCT) analogous to the Fast Fourier Transform (FFT). Use of 
the FCT provides a simple and powerful formalism for detection of signals with variable frequency 
just as Fourier transform techniques provide a formalism for the detection of signals of constant 
frequency. In particular, use of the FCT can alleviate the requirement of generating complicated 
families of filter functions typically required in the conventional matched filtering process. We briefly 
discuss the application of the FCT to several signal detection problems of current interest. 
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I. INTRODUCTION 



The detection of periodic signals is a well-developed art. In contrast, the detection of signals with variable frequency 
is an active area of research in signal processing. Considerable progress has been made in recent years using a variety of 
time-frequency techniques which include wavelets, bilinear transforms, and Short Time Fourier Transforms (STFTs) 

Si- 

In this paper, we consider the detection of deterministic signals with unknown parameters. The case of deter- 
ministic signals with unknown amplitude, phase, frequency, and arrival time has been treated in the literature 
In this paper, we generalize to an arbitrary number of parameters and consider signals with a deterministic, but 
parameterized, frequency evolution. We call these "variable frequency signals" . Specifically, we consider signals of 



the form: 



h (t) = i A (-t)cos(4>(t, A , A M )) < t < T 
s ^ ' otherwise ^ ' 



where (p(t, Ao, Am) is real and the {Ao,...,Am} represent various parameters which describe the phase evolution 
(e.g. frequency, frequency derivative, etc.). In general, </>() may depend non-linearly on time. For this paper, the 
"instantaneous frequency" of h s must be a well-defined quantity and the frequency evolution of the signal must be 
well-resolved in the data. Variable frequency signals include the class of signals usually known as chirps, i.e. signals 
which have a monotonically increasing or decreasing instantaneous frequency. Such signals appear in many contexts, 
such as almost-periodic signals with a small frequency drift or periodic signals emitted from accelerated systems. 
Chirps are discussed often in the literature of signal processing |^J4|,^] • We will concentrate on chirp signals in order 
to simplify the description of the chirp transform algorithm. 

A standard technique the for detection of signals in the presence of noise is the "matched filter" technique ||,|| . 
Detection of a signal h s (t) in the presence of white noise in a data stream h(t) of length T is based on the matched 
filter output: 

S = { dth s (t)h(t) = [ air h fm (r)h(T - t) (2) 



where hfut(t) is an optimal filter function in the time domain. For white noise, the optimum filter has an impulse 
response given by hfn t (r) — h s (T — t), in the interval < r < T. To detect a signal beginning at time to in a data 
stream of arbitrary length we compute the quantity: 

/>oo 

5(t ) = / dr h fm (T)h((t +T)-T) (3) 
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2 



where h(t) is zero outside a finite region of interest. When colored noise is present, it is conventional to work in 
the frequency domain. For matched filtering of real signals with an unknown arrival time, one constructs a signal 
estimator of the form 



S(t ) = 4i?e 



df 



Sh(f) 



where h(f) is the Fourier transform of the signal plus noise, h(t), defined as 

h(f) = / dt h(ty 2 ^ , 



(4) 



(5) 



h s (f) is the Fourier transform of the signal waveform and Sh{f) is the one-sided noise power spectral density. An 
optimum filter output is calculated for each realization of h s from Eq. [TJ using different values of the parameters. 
The approach of matched filtering thus requires the construction of a "dense" set of signal waveforms which cover the 
parameter space of possible signals. 

Other approaches to the variable frequency detection problem include techniques based on either synthesizing a 
multichannel filterbank Q or resampling the data at a variable rate These techniques are used widely in the 

radio pulsar community and, like the conventional matched filtering approach described above, differ significantly 
from the algorithm presented here. 

Just as the Fourier transform can be considered a form of matched filtering using a dense set of sine and cosine 
functions appropriate for periodic signals of unknown amplitude and phase, we wish to define a transform which 
performs matched filtering for a dense set of variable frequency signals with unknown parameters. By using the 
Fourier transform for periodic signals, we avoid the need of explicitly storing and computing the individual filter 
functions, i.e. the sine and cosine functions in the time domain. Analogously, the appropriate transform for variable 
frequency signals will avoid the need of generating a large set of filter functions and will provide a prescription for 
densely covering the set of possible signal waveforms. We will informally call the transform for variable frequency 
waveforms a "chirp transform" . 

The term "chirp transform" has been used elsewhere in the literature. For instance, Oppenheim et al. JToj] describe 
a "chirp transform algorithm" (CTA) which is a special case of the "chirp-z transform" (CZT). The chirp-z tranform 
is well-known and can be used to evaluate quadratic chirps. The method described in our paper is general and not 
constrained to quadratic chirp functions. We call the algorithm described in this paper the "fast chirp transform" 
(FCT). 

The techniques discussed in this paper appear to be related to filter bank design, to wavelet analysis, and to STFTs. 
In fact, the chirp transform can be viewed as a prescription for coherently adding the outputs of a bank of variable- 
length STFTs with a particular time-domain relationship. Thus, the fast chirp transform may already exist in some 
other formalism in the signal processing literature, but we are currently unaware of it. The recent work of Schutz ]29| ] 
and Williams and Schutz [fll| describes an approach which in similar in several aspects, but differs in that constant 
length STFTs are used rather than variable-length transforms. 

In section ||, we give the discrete forms of the matched filter outputs and discuss how these may be expressed 
in the form of the discrete analog to generalized Fourier integrals. In section III we derive the Two-Parameter Fast 



Chirp Transform that can be used to evaluate the discrete matched filter expressions. In section IV we generalize the 
definition of the FCT to an arbitrary number of parameters. In section ^ we briefly discuss several applications of 
the Fast Chirp Transform. 



II. DISCRETE MATCHED FILTERING 

The discrete forms of the time-domain and frequency domain matched filter outputs (Eq. || and ^) are given by: 



No-l 



<S(A , A M ) = ^ h s (j,X ,...,X M )h(j) , and 

3=0 



(6) 



S(jo,Ao,...,Am) = -rrRe 
i\ 



Hk )K(ko,\o,---,XM)e-^ k °/ N ° 



E 

k =0 



Sh(k ) 



(7) 
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Appendix [A| shows that for variable frequency waveforms, h s (t) = A(t) cos((/>(i, Ao, Am)), the discrete matched filter 
outputs can be expressed as generalized Fourier integrals, and in discrete notation take on a particularly simple form: 



<S(Ao,...,Am) = Re 



N -l 

Eft 

io=o 



(8) 



i\ 



JVo-l 

E '"■ 

,fco=0 



(9) 



where Qj and 7Yfc can be considered as the tim e or frequency series to be transformed, and $(fco, joi Ao, Am) is a 



real phase function of the form given in Eq. A9 



In the next section we define the Chirp Transform and show how it can be used to evaluate discrete transforms of 
the type shown in Eq. and |. Appropriate forms of the FCT will replace both the forward and inverse transforms 
contained in equations gfand |E . An inverse FCT is not applicable to the detection problem and will not be considered in 
this paper. The above formulation of the matched filtering process included the starting phase as a search parameter. 
An alternative approach is to convolve the signal with both the in-phase, cos(0), and quadrature-phase, sin(0), filters 
and then sum the squares of the results. This formulation is independent of the starting phase. In order to make the 
chirp transform applicable to both formulations and to a wider class of problems in general, the Re[] operator will 
not be included in the definition. If necessary, this operator can be applied after the transform. 



III. THE 2-PARAMETER FAST CHIRP TRANSFORM 
A. Example: The Quadratic Chirp (Linear Frequency Drift) 

As an initial example, we consider the problem of the detection of a quadratic chirp, e.g. a signal of the form 
h s (t) — A(t) cos(27r(/t + ^fi 2 )). We wish to detect this signal by matched filtering with a dense set of quadratic 
chirp waveforms. This requires evaluation of sums of the form in Eq. §. Note that if / were zero, the signal waveform 
would be periodic and we might discretely sample the input signal and then compute a power spectral estimate using 
the Fast Fourier Transform (FFT). 

Here we define a Fast Chirp Transform (FCT) for the quadratic chirp analogous to the FFT. This definition will 
be generalized to an arbitrary parameter frequency waveform in the next section. For simplicity we first define the 
Discrete Chirp Transform (DCT) for the quadratic chirp in analogy with the Discrete Fourier Transform (DFT): 

No-l 

H{koM} = E ^ e l2 " (fco(3o/Wo)+fcl0o/Wo)2) (10) 

Jo=0 

where hj is the discretely sampled data. The quadratic nature of the chirp is specified by the term (jo/No) 2 in the 
exponential. If k\ is zero, we have the usual DFT. 

To derive the Fast Chirp Transform (FCT), we begin by breaking the interval {jo = 0, ...,Nq — 1} into a set of 
contiguous sub-intervals 

iVi-i j™ in (,ji+i)-i 

01=0 jo=j™ in (ji) 

where jo nm (jx) * s the lower boundary of each of the intervals. The requirement that the intervals be contiguous implies 
io" n (j'i + 1) = jh rim (h) + n o(ji)j where no(ji) is the number of points in the interval. 

The boundaries are chosen as follows: we demand that the term ki(jo/N ) 2 change by no more than tt for values of 
k\ appropriate for the quadratic chirp signals being considered. If the term ki(jo / Nq) 2 changes by more than tt over 
a single sample, then the signal is not considered to be finely enough sampled to resolve the phase evolution. This 
limitation is analogous to the Nyquist sampling limit for periodic signals. Requiring that the term ki(jo/No) remain 
relatively constant over a sub-interval allows us to approximate the DCT as: 
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H 



{fc ,fcl} 



Ni-1 



Jo 



l Ui)/Nof) 



30- 



E 



'Jo' 



J2TTko(jo/N ) 



(12) 



and we note that the second summation can now be computed as a Fast Fourier Transform (FFT). We further demand 
that the term fei(j™"(ji)/iVo) increment by a constant amount from one sub- interval to the next. The requirement 
that the increment be less than 1/2 for the maximum value of fci, k™ ax , specifies Ni = 2k" lax , with the "Nyquist" 
restriction that k™ ax < Nq/4 for the case of a quadratic chirp. We can then write: 



H 



JVi-i Jo™* n (ji+i)-i 

J'i=° JO =j™ in Ui) 



(13) 



where we have anticipated that we will evaluate the sum for integral values of fci using a standard FFT. Finally, 
we express the inner sum in standard FFT form by extracting a phase factor and we define the 2 parameter FCT, 

C{fc ,fci}7 as: 



E 

J'i=0 



J2irko(j™ 



no(ji)- 

'(ji)M)) 

30=0 



jo+io" n (ji) 



(14) 



As before, no(ji) is the number of points in the interval. We call fco and fci the "conjugate variables" of the linear and 
quadratic terms, respectively, in the same sense that frequency and time are a conjugate variable pair. Note that both 
summations can be implemented using FFTs and thus we are justified in calling the transform a fast chirp transform. 
The non-linearity of the quadratic chirp is absorbed into the specification of the boundaries of the contiguous intervals. 
This is the key concept of the FCT. Note that the definition of CW^j is general and does not depend explicitly 
on the quadratic nature of the phase evolution. The same definition will apply to other non-linear phase evolution 
functions. 

We may also write: 



C 



{fe 0l fci} 



JVi-1 

E ei 

31=0 



!7rfei(ji/ATi) r 

^{feo 



(15) 



where we have used the notation Cf kn.h \ to indicate that it is a partial transform, i.e. transformed over one index, 
jo, but not the other, j\. Equation Hq illustrates an interesting feature of the FCT, which is important in imple- 
mentation considerations. Although the number of points in an interval, rio(ji), may be small, the partial transform, 
Cj^ jj, requires evaluation at a large number of values of fco, for example at iVo values of fco- This is equivalent to 
calculating the oversampled FFT of the individual intervals with an oversampling factor of iVb/no (ji ) . The problem 
of computing the C^k ,ji} thus reduces to the problem of estimating the oversampled spectrum of each interval. The 
oversampled spectrum may be calculated exactly using the Fractional Fourier Transform (FRFT) . Fast methods have 
been developed for evaluating FRFTs jl^] and it can be shown that the calculation of the values of the individual 
oversampled FFTs that enter into C{k ,ji} m Eqs. [l4| and |l5| requires 0(A r olog2no(ji)) operations to be computed 
exactly. The computation required to calculate the entire set of FRFTs can be shown to be 0(NiNolog2(No/Ni)). 
The phase factors, exp (i2nko(j^ lm (ji) / Nq)) in Eq. [l4| are easily computed as part of the same formalism. 

Taking into account the evaluation of the inner and outer sums separately, the number of operations required for 
the evaluation of the FCT can be shown to be: 



AV < O(AWog 2 A ) 



(16) 



which is at least as efficient as the matched filter approach. The inequality indicates that the evaluation requires less 
computations if approximations are used in evaluating the oversampled FFTs or if a coarser sampling of the FCT 
over fc values is used. We will discuss the issue of computational efficiency further in Section III D below. 



B. Accuracy of the Approximations 



How accurately can the FCT approximate the discrete matched filter output? There are three types of approxima- 
tions to be considered: (1) Possible use of the Stationary Phase Approximation in deriving a form for the matched 



5 



filter output; (2) The constant phase approximation used to derive the discrete form of the FCT; and (3) Possible 
use of approximations in calculating the oversampled FFTs that enter into the FCT. We will consider each of these 
approximations in turn. 

(1) There may be some error due to the Stationary Phase Approximation itself (Appendix ^) if this is used to 
approximate the Fourier transform of the signal waveform. This is not inherently part of the FCT, and we will not 
discuss this approximation here since it depends on the specific waveform being analyzed. However, we note that the 
Stationary Phase Approximation can be extremely accurate in practice [ OLjijJ . 

(2) The error due to the constant phase approximation (Eq. [l2|) in the FCT itself is obviously dependent on the value 
of the conjugate variables at which the FCT is evaluated. If the value of k\ is small, the value of the increase in the 
quadratic phase term as j\ goes to j\ + 1 is also small, and the approximation will be very accurate. By analogy 
to the Fourier transform, we expect the "frequency response" of the output of the matched filter computed by the 
FCT to behave similarly to the frequency response of a power spectrum computed with an FFT. Specifically, just as 
there is a roll-off in power for some periodic signals near the Nyquist frequency, there will be a roll-off in the accuracy 
of the FCT approximation to the matched filter output for signals near the Nyquist limit of the conjugate variable, 
namely fci (Nyquist) = N\/2. For values of k\ well below Nyquist, we expect the accuracy of the FCT approximation 
to the matched filter output to be very good. Furthermore, if we wish higher accuracy in the FCT approximation, 
we can employ the same techniques used in Fourier analysis, namely, we can "oversample" the FCT. This can be 
accomplished by "zero-padding" of the outer (ji) FFT. 

(3) Finally, there may be errors due to the possible use of approximations in calculating the oversampled FFTs. As 
indicated in Eq. [l6], the exact calculation of the oversampled FFTs requires O(NiN \og2N ) calculations. However, 
this assumes that each interval is oversampled by a factor A /no(ji)i which can be quite large. Typically, oversampling 
factors of 2 2 — 2 3 are sufficient for an accurate approximation. If N± is large, the computational requirements may 
thus be reduced considerably by using such approximations (i.e. 0(Aolog2-/Vo) rather than 0(Ai Aolog2 No)). A 
quantitative discussion of oversampling approximations is beyond the scope of this paper. 

We also remark that the FCT formalism lends itself naturally to a variety of hierarchical search approaches. For 



instance, consider Eqs. |1J and 15. The outer sum is a coherent addition of the contributions from the individual 
intervals. In order to implement a fast hierarchical search, we could simply take the magnitudes of each C{k ,ji} an d 
add them as an incoherent sum over ji, giving a measure of the incoherent power as a function of fco. Values of ko 
with significant incoherent power could then be examined in more detail by performing the coherent summation over 
ji. This could lead to a dramatic decrease in the required number of computations, depending on the threshold set 
for the incoherent power summation step. 



C. Implementation Considerations 

As an initial trial, we have implemented the two-parameter FCT and tested it on se veral types of waveforms, 



including the quadratic chirp discussed above and the "Newtonian chirp" discussed in Sec. V A below. We have used 
two implementations, one which uses a fixed length oversampling of the inner FFT, and one which uses a pre-packaged 
2D FFT algorithm. 

In this section we will describe in more detail the 2D FFT implementation in order to give further insight into the 
details of the 2-parameter FCT algorithm. The 2D FFT implementation is extremely simple to code although it is 
not the most computationally efficient. It will be shown in the next section that even so, it is nearly as efficient as 
the brute-force matched filtering method. 

The 2D FFT implementation of the FCT packs the initial one dimensional data array, h.j Q , into a sparse two 

dimensional array hj j 1 . The packing will be determined by the j™ m (jx) array which is ultimately defined by the 
phase function 4>(jo). Once the data are packed appropriately, the FCT is calculated using any pre-packaged 2D FFT 
routine. For the rest of this discussion, it is assumed that Nq and A^i are compatible with the 2D FFT routine being 
employed. This usually means that these lengths are a power of two. 
From equation |l3|, the two parameter FCT may be written as 

Ni-l j^ in (ji + ^)-l 

CKM = J2 e l2 * klUl/Nl) J2 ^oe l2 " MWJVo) . (17) 



By defining a two dimensional Nq x Ni array such that 
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t _ J K j mm (ji) < io < ir n (i! + 1) n ^ 

J0jl ~ \ otherwise ' 1 J 

equation ^ becomes 

iVl-l JVo-l 

-EE /^e^W^+W*.). (19) 

31=0 3*o=0 

As promised, C{fc ,k ± } is the 2D FFT of the sparsely packed array hj j 1 . For the case of a monotonic phase function 
that satisfies (f)(0) — and 4>(Nq) — Ni, the array boundaries are given by 

io min O'i) = r 1 (ii)- (20) 

Note that we have the freedom to rescale <j>(j) so that N\ may be chosen arbitrarily. 

The number of operations needed by the FCT implementation using a 2D FFT is the same as that needed by a 
two dimensional FFT of order Ao xJVi: 

N fct ~2D = O (NqNx log 2 N Nx) . (21) 



D. Computational Efficiency 

For a matched filter operation using individual filters, there is one Fourier transform to perform for each filter, thus 
the number of computations, N m f, is of order: 

N mf = O (JVxiVolog2iVo) (22) 
where Ai is the number of filter functions needed to cover the space of possib le waveforms, and Nq is the number 



of samples in the time series or frequency spectrum. As shown in section III A using the fractional fouricr transform 
(FRFT), the FCT computation can be of order 

A fct—FRFT < O (Ax A log 2 N ) . (23) 
Two approaches offer po tentia l computational gains: (1) Reducing the oversampling factor in computing the inner 



FFT (discussed in Section HI A ), and (2) Relaxing the requirement of full resolution in the fco variable. Reducing 
the oversampling factor can change th e inn er sum computation from an O (AoAi log 2 (Ao/Ai)) calculation to an 
O (Nq log 2 Nq) calculation (see Section III A| ) . If N% >> 1 the total calculation of the FCT will then be of order 



O (AoAilog 2 Ai). This is potentially a factor of O (log 2 Ao/log 2 N±) more efficient than the conventional matched 
filtering technique. We emphasize that the detailed coefficients in front of the scalings are not yet known for compu- 
tationally efficient implementations. 

Significant computational gains are also potentially available by relaxing the requirement of full resolution in the ko 
variable. Not all problems require the high ko resolution that the previously discussed FCT implementations deliver. 
While the power is very localized in the ko variable when the value of k% is that of the actual signal(and vice- versa), the 
power can be significant for values of ko and k\ which simultaneously deviate from the actual signal values. The reason 
for this is that the deviation in ko can be compensated for by a deviation in k± , providing a reasonable correlation of 
the matched filter template with the actual signal. The determination of the optimum sampling resolution is closely 
related to calculation of the ambiguity function |j| . Discussions of techniques for calculating the ambiguity function 
have been given by Owen ]Tq] , Mohanty and Dhurandhar and Owen and Sathyaprakash |l7|] for special cases of 
the types of chirp functions considered here. 



The FCT may be evaluated at lower resolution by reducing the order of the 2D array hj j t (Sec. HID) from 
Ao x Ai to (eAo) x Ai where e < 1. As long as (cAq) x N\ > No, we can still pack the original Ao data points into 
this smaller array. The resulting FCT will have a coarser fco resolution but it will take 0(eAiAolog 2 Ao) operations 
to perform (again neglecting terms of order log 2 iVi). Alternatively, in the straightforward evaluation of the FCT 
contained in Eq. the sums over ji are simply evaluated at a decimated set of k Q . 

We remark that sampling the FCT at lower resolution can be used as part of a hierarchical algorithm, i.e. the FCT 
is first calculated at lower resolution and values of fco are identified having excess signal strength. The FCT is then 
evaluated on a finer grid of fco values near those values of fco exhibiting the excess signal. 
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E. Selection of the Range for Evaluation of the FCT 



We can select ranges of ko and fci to be evaluated by a process analogous to a heterodyne operation on a periodic 
signal. Selection of the range of ko is done by the usual process of down-conversion and low-pass filtering and we will 
not consider it further here. For the parameter k\ the process is different in that no low-pass filtering is required, and 
because the range of k\ can depend on ko, i.e. the conjugate parameter ranges need not be independent. This provides 
considerable flexibility in determining the shape and volume of parameter space that can be efficiently searched. 

Suppose we wish to compute the FCT for a range of k\ centered on k™ ld . Re- writing the 2-Parameter FCT including 
the term, we obtain: 



JVo-l 



Mr, , , - V h ■ p i^(koUo/No)+(ki+kr d (ko))(jo/N ) 2 ) 

-"{feo.fci} — 2-~i 3a 

3o=0 



(24) 



where we have explicitly allowed k™ ld to depend on fc - As before, assuming that k\ + k™ ld is sufficiently small, we 
take the quadratic phase term out of the inner summation to obtain: 



ATi-l 



H{k a M} ~ C {k M) - E 



a i2irfci(ji/JVi) 



ji=0 



«oO'i)-l 

^(koU^ij^/^+kr^koKh/N,)) h . ( 

30=0 



i27rk (j /N ) 



To further simplify the notation for the outer sum, we define: 

e(ko,h,jr n ,kT id ) = k (jr n (ji)/No) + fcr d (fco)(jiM) 

to obtain 



JVi-l 

c {k0 , kl } = E ^ {h/Ni) 

31=0 



"o(ji)-l 

J2ne(k ,j™ in ,k? id ) h m . n i2nk Uo/N ) 

Jo=0 



(25) 



(26) 



(27) 



ft is important to note that N\ can be chosen to be any integer less than or equal to the "Nyquist limit", No/2. This 
is due to the fact that we can choose the boundaries, j™ m (ji), of the intervals arbitrarily as long as the phase change 
from the quadratic term is kept sufficiently small. Thus, we are free to choose any range of k\ around k™ ld as long as 
the size of the range does not exceed the Nyquist limit. 



IV. GENERALIZED DEFINITION OF THE FCT 



Because the non-linearities of the phase evolution enter into the FCT only in the boundaries of the sub-intervals, 
we may generalize the definition of the FCT to the case of a sum of non-linear phase evolution terms. In particular, 
consider the DCT: 



JV -1 

#{fc ,..,fc M -i} = E ^'o ex P 
3a=0 



(28) 



where the {(f> p (jo) ■ p — 1, M — 1} are a set of parameterless, non-linear phase functions. The corresponding FCT 
is: 



C 



{k a ,...,k M -l} 



JVm-i-1 

E e 

j M -i=0 



(29) 



x e 



j™'"( J!r ..j M _ 1 j 

i2irk\ ( — ^ ) 



x e 



, , ■ ■3M-l) . 
*2vrfe ( jv^ ) 



™1 (j2vjM-l)-l 

E < 

31=0 
n 0'i,...,jiif-i)-l 



,*27rfei(-^-) 



^ e ^ko(^ >h 



jo+jy in (ji,...jM-i) 



io=o 
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where the j™ m {j p+ i, -3m-i) and N p (ji, ...Jm-i) are specified by the phase functions {0 p (jo) ■ P = h —,M - 1}, 
and the {N p : p = 1, M — 1} are determined by the maximum allowed values of the {k p : p = 1, M — 1}. As 
before, the intervals specified by the j™™(ji, ■■■JM— l) are contiguous. Thus for each non-linear phase function 4> p (jo) 
there is a corresponding "conj ug ate variable", k p . 



In analogy with equations p_8| and 19, the M parameter FCT may be written in the form of an M-dimcnsional 



discrete Fourier transform by defining the matrix hj 0t .,_ t j M _ t such that 



h j0 ,..., jM ,. = | > ^ jpT < * + 1 fOT a11 P e t 1 ' M " !] . (30) 



otherwise 
With this definition, equation ^9] becomes: 



Nm-i-1 N a -1 



Note that the interval boundaries, j™ m , may be determined from equation [30|. 

Finally, we consider how the range for evaluation of the generalized FCT can be specified, analagous to Eq. |2?] 
for the 2-parameter case. To specify the region for which the FCT is to be evaluated, we add a term to each of the 
exponential terms in Eq. As for the 2-parameter case, equation we define: 

Q P (k P -i, j P , j™?, k™ d ) = Vi(j P m -T/Wp-i) + k™ id {j p /N p ) (32) 

where j™" is a function of {j p , ...,j'm-i} and fc™ d can depend on {fco,fci + k™ ld 1 ...,k p -i + fc™i}- The parameter 
space searched by the FCT will then be k p nld ± N p /2. This provides considerable flexibility in determining the shape 
and volume of parameter space that can be efficiently searched. We can now write: 

IVm-1-1 , 



C {k0 ,..., kM -,} = 2J e (33) 

»M-2(jM-l) — 1 



3m-2=0 



ni(j2,-"jJM-i)-l 

x el 2 7r e 2 (fe 1J 2jr" , ^? i * d ) e i2 ^fei(^-) 

ii=o 

naO'i,— Jm-i)- 1 
l 27re 1 (feoJiJo ml "^r" d ) V ^ ko[ ^ ] h 

A C ^ "•j0+J™ , "0l:---jAf-l)- 

J0=0 



(34) 



V. DISCUSSION: APPLICATION TO DETECTION OF VARIABLE FREQUENCY SIGNALS 
A. Detection of Gravitational Waves from the Binary Inspiral of Neutron Stars and Black Holes 

One of the primary goals of the new generation of laser interferometric gravitational wave detectors is the detection 
of gravitational waves from the binary inspiral of compact objects, specifically neutron stars (NSs) and black holes 
(BHs). There is a large literature written on the subject of matched filtering for detection of gravitational waves 
using laser interferometers (see e.g. |L8| p0[ ) . The matched filtering techniques are based on Eq. ^ where h(f) is the 
Fourier transform of the gravitational strain generated from the differential output of the interferometer, the h s (/) 
are the Fourier transforms of theoretically generated binary inspiral signal waveforms, and Sh(f) is the measured 
power spectral density of the interferometer. 

A significant amount of work has gone into the calculation of binary inspiral waveforms (called "templates" ) , the 
spacing of such templates to achieve near-optimal sensitivity, and the cost of generating such templates in terms of 
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compute cycles and storage requirements pl| , [l5| , p^7| , p2| , p3t . Current matched filter techniques require thousands to 
tens of thousands of templates to cover the space of expected waveforms depending on the mass range of the binary 
systems considered. 

The method of chirp transforms described here does away with the requirement of generating thousands of individual 
templates and provides a natural way to cover the space of allowed waveforms completely. To apply the chirp transform 
to the binary inspiral problem, we make use of the stationary phase formalism. Droz et al. |13|] have shown that the 
stationary phase formalism can be used to provide an accurate approximation to the Fourier transform of the time- 
domain waveforms of inspiraling binaries as calculated in the "Newtonian" approximation. This is essentially an 
application of the stationary phase approximation discussed in Appendix [X] to the case of gravitational waveforms. 
Damour et al. (24||l4|] have shown that the binary inspiral waveforms can be accurately calculated using the SPA and 
an alternative formalism based on "P-approximants" . They note that care must be taken in the treatment of the 
termination of the waveform at the time of the final plunge and merger. 

In order to illustrate the use of the FCT in gravitational wave detection, we discuss the example of a "post- 
Newtonian" (PN) expansion. The canonical PN stationary phase waveforms for binary inspiral up to 2PN order are 
of the form @: 

*■<"- (^) 1 ' 2 (^) 1,3 ^ 7,6T - 1, ^i i *( f »i < 35 ' 



*(/) - 27r/t c - 2<f> c - tt/4 



+ 



128ry 



x -5/3 + 



/3715 55 



V 756 



/ 15293365 27145 



V 508032 



504 



-v - 



3085 



72 



-V X 



- 16?rx" 2/3 

-1/3 



where, for simplicity, we consider only one polarization. The variables have been defined as usual (M is the mass 
of the sun, T© is GM© /c 3 and has a value of approximately 4.925 x 10 _6 sec, Mtot is the total mass of the binary 
system, 77 = /i/M tot , [i is the reduced mass of the binary, t c is the time of coalescence, 4> c is the phase at coalescence), 
and we have defined 

X = - 1 g_0/ (36) 

It can be seen that the stationary phase waveforms have frequency dependent amplitudes and phase functions that 
are expansions in powers of the frequency, /. In particular, 

*(/) = a + 2ntJ + A /- 5 / 3 + X 1 f~ 1 + A^r 2/3 + A^ 1 / 3 . (37) 

where a is a phase constant and the X x are coefficients of the frequency expansion which depend on Mtot, M©, T©, 
and 77. 

In order to apply the FCT, we construct the discrete version of the matched filter output, Eq. 0, where h(ko) is the 
Fourier transform of the discretely sampled interferometer strain output, h s (ko) are the stationary phase waveforms 
given in Eq. [?5| above, and Sh{ko) is the noise power spectral density of the interferometer. 

The FCT is then used to evaluate the matched filter, with a resulting output, 

C{t c: Ao,A 1 ,A 1 . 5 ,A 2 } (38) 

In this expression, 2irt c is the conjugate variable of the linear / term; Ao is the conjugate variable to the Newtonian 
term, / _5 ^ 3 , Ai is the conjugate variable to the 1PN term, etc. Numerous considerations arise in selecting 

the search ranges for the conjugate parameters. As discussed by several authors |2^Jl7[] , spinless low-mass binaries 
should be reasonably well detected by a three parameter search over {t c , Ao, Ai}. The conjugate parameters Ao and 
Ai fulfill a function similar to the parameters r and t\ that appear in the literature (e.g. J22|;[l7]]) which represent the 
Newtonian and 1PN contributions to the time to coalescence, respectively. Thus, in the case of spin-less low-mass 
binaries, A^ x 5 2 } can be considered functions of A{o,i}- Owen and Sathyaprakash [jlTj point out that T{ ,i.5} may 
be more convenient search parameters. Hence, A{ 12 } become functions of A{ .i.5}- In this case, the "heterodyne" 
approach (see the discussion following Eq. |3l|) can be used with the dependent parameters to reduce the search space 
to that of spin- less 2PN waveforms. 
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To search for binaries with spin, additional independent parameters will be needed and thus it will be useful to 
search in a range around the spin- less PN expansion (or other expansion). This can also be accomplished using the 
method described in the discussion following Eq. |l[ The technique will be particularly useful for massive binaries 
for which spin interactions could be significant. An important step will be to estimate limits to the range of the 
conjugate variables in the FCT analysis due to spin effects. The FCT then provides a formalism for searching the 
complete waveform space, even if the exact waveforms are not known. 

It will also be quite useful to enlarge the search region beyond the space physically accessible by astrophysical 
binary systems. While no binaries are expected outside the physically accessible regions, it is important to study the 
characteristics of noise signals in regions close to the physically accessible regions. The FCT formalism provides a 
straightforward way to tailor the analysis to a range of search regions. This, of course, is also possible with conventional 
template-based techniques. 

The FCT formalism may be useful for expansions other than post-Newtonian. In particular, we are very interested 
to see whether the FCT formalism can be applied to the P-approximants discussed very recently by Damour et 
al. 01. Also, as we remarked earlier, the FCT lends itself naturally to hierarchical approaches for binary inspiral 
detection. We note in particular the recent paper by Tanaka and Tagoshi which discusses efficient hierarchical 
search algorithms which have several similarities to the general FCT algorithm. 

In summary, the use of the FCT for detection of the chirps from gravitational waves has several attractive features. 
First, no explicit calculation and storage of gravitational waveforms is required for the analysis. Only the order of 
the PN expansion, the power-law exponents appearing in the expansion, and the range of the search parameters is 
important. Second, waveforms with perturbations on the phase evolution such as those due to spins can be detected 
even if the exact waveforms are not known since the FCT can be used to search completely an arbitrary region of 
parameter space. The only requirement is that the perturbation not involve significant terms beyond those in the 
expansion considered for the FCT. The FCT may therefore be very useful in the search for BH-BH coalescence where 
the waveforms are not precisely known p5[ , or for sources to be detected by the space- based gravitational detector, 
LISA, where the waveforms are also only approximately known and phase perturbations are likely to be present. 
Finally, the FCT formalism can be used to investigate the characteristics of noise signals in the neighborhood of 
expected signals from binary inspiral. 



B. Detection of Rotating Neutron Stars 



1. Acceleration Searches 



Pulsars are rotating neutron stars that spin with periods in the range of approximately 1 ms up to hundreds of 
seconds. Pulsars are detected primarily at radio and x/7-ray wavelengths. In the future, rotating neutron stars may 
also be detectable as sources of gravitational waves. Detection of pulsars usually employs Fourier transform techniques 
to find the periodic pulses. However, several effects complicate the detection of pulsars and cause the pulsations to 
deviate from being strictly periodic. For instance, the emission from pulsars in compact binary systems is Doppler 
shifted causing a frequency variation on the time scale of the orbital period. Likewise, the earth's rotation and orbit 
can induce frequency and phase variations that are dependent on the position of the source on the sky. Rotating 
neutron stars can also have non-negligible spin down effects, especially if the neutron star is young. Any of these 
effects can be important at both radio and x/7-ray wavelengths depending on the length of the observation. They 
are also likely to be important in future searches for gravitational wave emission from rotating neutron stars due to 
r-modes, or from older rotating neutron stars because of the earth's orbit and rotation. 

In the past, so-called "acceleration searches" have been employed to detect pulsars with slowly varying frequency [||. 
These are essentially matched filter techniques implemented either with templates, or equivalently, with "stretching" 
of the time or frequency variable. This requires individual matched filter operation, one for each discrete acceleration 



trial. The FCT analog is that of the quadratic chirp analysis discussed in Sec. [II. The FCT also provides a natural 
extension to searches beyond quadratic (acceleration) effects. 



2. Dispersion Measure Searches 



Radio radiation emitted by pulsars travels through a diffuse interstellar plasma known as the interstellar medium 
(ISM) before reaching detectors on Earth. The dispersive properties of the ISM cause individual radio pulses to 
broaden in time. This dispersion broadening will reduce the chances of detecting a given pulsar signal. The magnitude 
of the dispersion effect is measured by a quantity called the dispersion measure (DM). If the DM is known, the 
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dispersion effect can be removed from the pulsar signal using standard digital signal processing (DSP) techniques. 
When searching for new pulsars, the DM is rarely known and systematic searches must be performed both in DM 
and in pulsar period. 

The effect of ISM dispersion may be removed from the received signal by applying the following transformation 
©1: 

/OO 
S r (f)e l27rft e l27TDMli > u) df (39) 
-oo 

where S r (f) is the Fourier transform of the received signal and S(DM, t) is known as the dedispersed signal. When 
searching for new pulsars, one must generate several time series corresponding to different values of the DM. The 
above equation shows that S(DM,t) is simply a chirp transform of S r (f). The FCT provides an efficient way to 
generate S(DM, t) for several values of DM. Each of these time series can then be searched for periodic signals. 

When searching for pulsar signals, one typically "detects" the total power in the signal by calculating P(DM,t) = 
\\S(DM, t)\\ 2 and then averages over a small window of time. Each time series is then searched separately using various 
pulsar detection techniques. The structure of the FCT points to the possibility of a slightly different technique. Rather 
than searching each time series separately, one f irst calculates P(t) = J2dm P(DM, t) and then searches this time 



series for possible pulsar signals (see section piB|) . Using the property that the sum of the squares is conserved under 
a Fourier transform, the second set of FFTs in the FCT does not need to be performed in order to calculate P(t). 
Thus, a highly efficient intermediary chirp transform can be used instead of the complete FCT. 



VI. SUMMARY 



We have described an algorithm for the detection of signals with variable frequency. Standard detection algorithms 
use matched filtering techniques which require both the computation of a large set of task specific filter functions 
and a prescription for densely covering the set of possible signal waveforms. The Fast Chirp Transform proposed in 
this paper automatically precludes the need to generate specific filter functions since standard FFTs can be used in 
the implementation and the FCT immediately provides the prescription for densely covering the waveform parameter 
space. 

The FCT for a two parameter chirp was defined and then generalized to N parameters with arbitrary phase functions. 
A straight forward implementation of the FCT was discussed and it was shown to be comparable in efficiency with 
the brute-force matched filtering approach. Several approaches to achieving even better computational efficiency were 
also discussed. 

The efficient detection of variable frequency signals has a large number of practical applications. Of considerable 
interest to the authors is the detection of gravitational waves from NS and BH binary systems and the detection 
of radio waves from pulsars. Another obvious area of application is radar-sonar signal processing where target or 
transmitter motion can cause Doppler frequency shifts in the received signal. Other potential areas of application 
include communications and image processing. A more detailed description of the FCT and its application to the 
above problems will be the subjects of future work. 

We acknowledge informative and pivotal discussions with Ben Owen and B. Sathyaprakash. We also acknowledge 
helpful and stimulating discussions with Alessandra Buonanno, Tibault Damour, Scott Hughes and Albert Lazzarini. 
This work was supported in part by the LIGO Laboratory under cooperative agreement NSF-PHY-92 10038 and by 
grant NSF-PHY-9970877. This paper is LIGO document LIGO-P000003-00-R. 



APPENDIX A: THE STATIONARY PHASE APPROXIMATION AND MATCHED FILTERING 



1. Frequency-Domain Matched Filtering 



We begin by showing how the Fourier transform of waveforms such as those of Eq. [j] can be approximated in 
a way that allows them to be expressed naturally in frequency-domain match ed filtering, Eq. ^. The Stationary 
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Phase Approximation (SPA) (see e.g. [g7J; see also Ref. |lq] and Refs. 
gravitational wave detection), provides a prescription for approximating the Fourier transform of a function of the 
form h s (t) = A(t) cos (4>(t)) (where A(t) and <fi(t) are real and 4>'(t) is positive): 



dt A(t) e 4/v, + (t) + / dt A(t) e-^-M 



(Al) 
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with i/j±(t) = 4>(t)/f ± 2Trt. If tf exists such that ip' + (tf) = or ip'_(tf) — 0. then tf is called a "stationary point". 
Considering positive / and positive 4>'{t), only the second integral in Eq. Al contains a stationary point. Hence, to 
leading order, we can write |2j: 



hs(f) 



dt A(t) e i(^/*-^(*)) 



(A2) 



for / > 0. Note that we compute h s (—f) using the Fourier transform property of real functions: h a (—f) = h*(f). 
If all derivatives of ip—{t) up to order p are zero at tf, then the Fourier transform of h s may be approximated by: 



with components given by somewhat complicated but straightforward expressions: 



Af) = A (t f ) 



pi 



i i/p 



f\¥ p) (tf)\ 



r(i/p) 



and 



(A3) 



(A4) 



*(/) - 2irft f -<j>(t f )±ir/2p 



(A5) 



where the sign of ir/2p is positive or negative depending on whether ip^f\tf) is positive or negative, respectively. In 
particular, for p = 2, the following approximation holds: 



h s (f) = J^AitfM'itfy-vwwt-Wf^w. 



(A6) 



The stationary phase approximation is accurate as long as the amplitude of h s does not vary too quickly compared 
to the time derivative of the phase, <p'(t), and the effect of the higher derivatives of <p(t) on the phase evolution are 
small compared to the effect of 4>'{t). 

Using the form given in Eq. A3, we can now rewrite Eq. 0. Gathering all the amplitude terms together, 



s h (f) 



we can express the matched filter output as: 

S(t ) = 4 Re / dfH(f)e-^ 
Jo 

where 

<*>(/) = *(/) - 2nft . 



(A7) 



(A8) 



(A9) 



Hence, the matched filtering operation in the frequency domain is expressed as an integral transform, specifically a 
so-called generalized Fourier integral. In analogy with the Discrete Fourier Transform (DFT), we can write this in 
discrete form as: 



S = —Re 

N 



N -l 
Lfco=0 



(A10) 



The summation can be computed as a FCT. 



2. Time-Domain Matched Filtering 

We note that for signal waveforms of the form, h s (t) — A(t) cos (<p(t)), the expression (Eq. ||) for time-domain 
matched filtering yields directly: 
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S = Re 



alt h(t)A(t)e zc f> w 



(All) 



Suc h sign als are of considerable interest and include periodic signals with frequency drift. The integral transform in 
Eq. All can be represented in discrete form in the usual way as: 



S = Re 



No-l 
jo=0 



where 

Gj = h(j )A(j Q ) 

and where h(ja), A(jo), and <p(Jo) are the discretely sampled values of the continuous functions. 
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